Journal  of  Power  Sources  192  (2009)  544-551 


ELSEVIER 


Contents  lists  available  at  ScienceDirect 

Journal  of  Power  Sources 

journal  homepage:  www.elsevier.com/locate/jpowsour 


pri 

SllbiidtS 


Assembly  pressure  and  membrane  swelling  in  PEM  fuel  cells 

Y.  Zhou,  G.  Lin,  AJ.  Shih,  SJ.  Hu* 

Department  of  Mechanical  Engineering,  The  University  of  Michigan,  Ann  Arbor,  MI  48109-2125,  USA 


ARTICLE  INFO 


ABSTRACT 


Article  history: 

Received  15  December  2008 
Received  in  revised  form  21  January  2009 
Accepted  28  January  2009 
Available  online  7  February  2009 


Keywords: 

Assembly  pressure 
Deformation 
Current  distribution 
Membrane  swelling 


Assembly  pressure  and  membrane  swelling  induced  by  elevated  temperature  and  humidity  cause  inho¬ 
mogeneous  compression  and  performance  variation  in  proton  exchange  membrane  (PEM)  fuel  cells.  This 
research  conducts  a  comprehensive  analysis  on  the  effects  of  assembly  pressure  and  operating  tem¬ 
perature  and  humidity  on  PEM  fuel  cell  stack  deformation,  contact  resistance,  overall  performance  and 
current  distribution  by  advancing  a  model  previously  developed  by  the  authors.  First,  a  finite  element 
model  (FEM)  model  is  developed  to  simulate  the  stack  deformation  when  assembly  pressure,  tempera¬ 
ture  and  humidity  fields  are  applied.  Then  a  multi-physics  simulation,  including  gas  flow  and  diffusion, 
proton  transport,  and  electron  transport  in  a  three-dimensional  cell,  is  conduced.  The  modeling  results 
reveal  that  elevated  temperature  and  humidity  enlarge  gas  diffusion  layer  (GDL)  and  membrane  inhomo¬ 
geneous  deformation,  increase  contact  pressure  and  reduce  contact  resistance  due  to  the  swelling  and 
material  property  change  of  the  GDL  and  membrane.  When  an  assembly  pressure  is  applied,  the  fuel  cell 
overall  performance  is  improved  by  increasing  temperature  and  humidity.  However,  significant  spatial 
variation  of  current  distribution  is  observed  at  elevated  temperature  and  humidity. 

©  2009  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Fuel  cells  are  promising  alternative  power  devices  due  to  their 
high  theoretical  efficiency  and  environmental  friendliness.  In  par¬ 
ticular,  proton  exchange  membrane  (PEM)  fuel  cells  are  attractive 
for  automotive  and  portable  applications  because  of  their  low  oper¬ 
ating  temperature,  fast  start-up,  and  low  emissions.  A  PEM  fuel  cell 
stack  consists  of  several  single  cells  connected  in  series  by  bipolar 
plates  (BPP)  which  provide  reactants  to  the  membrane  electrode 
assembly  (MEA).  Assembly  pressure  can  increase  the  overall  elec¬ 
trical  conductivity  of  the  gas  diffusion  layer  (GDL),  improve  contact 
resistance,  and  hence,  reduce  the  electrical  resistance  losses  inside 
a  cell.  Assembly  pressure  also  determines  the  contact  status  and 
stack  deformation  especially  that  of  the  GDL,  which  is  the  most 
deformable  component  in  a  PEM  fuel  cell  stack.  Under  the  land  of  a 
bipolar  plate,  the  GDL  is  compressed  by  the  assembly  force.  Under 
the  channel  area,  the  GDL  protrudes  into  the  channel  cavities.  The 
thickness  and  porosity  of  the  GDL  are  affected  under  compression; 
consequently,  the  mass,  heat,  and  charge  transfer  properties  are 
changed. 

It  is  well  recognized  that  the  GDL  can  influence  PEM  fuel  cell 
performance  significantly  [1].  However,  most  of  PEM  fuel  cell 
performance  models  do  not  consider  this  GDL  inhomogeneous 
compression  and  only  limited  research  has  been  conducted  to 
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address  this  issue.  Zhou  et  al.  [2]  developed  a  multi-physics  model 
to  investigate  the  effects  of  assembly  pressure  on  PEM  fuel  cell  per¬ 
formance  by  considering  contact  resistance  and  flow  resistance.  Sun 
et  al.  [3]  assumed  a  GDL  compression  ratio  and  analyzed  the  influ¬ 
ence  of  performance  and  current  density  distribution.  Hottinen  et 
al.  [4]  conducted  a  study  on  the  effect  of  inhomogeneous  compres¬ 
sion  of  GDL  on  the  mass  and  charge  transfer  in  PEM  fuel  cells  using 
experimentally  obtained  GDL  parameters  as  a  function  of  compres¬ 
sion  thickness.  However,  the  room  temperature  and  dry  conditions 
were  assumed  in  both  modeling  and  experimental  investigation  on 
GDL  deformation. 

Room  temperature  and  dry  conditions  do  not  reflect  the  real 
PEM  fuel  cell  operating  conditions  since  most  fuel  cells  operate  at 
elevated  temperatures  and  100%  relative  humidity  (RH).  Elevated 
temperature  and  high  RH  influence  PEM  fuel  cell  polarization  losses 
in  many  ways  including  catalyst  activity  [5],  membrane  mechan¬ 
ical  and  electrical  properties  [6],  and  gas  transport  [7],  etc.  In 
addition,  GDL  and  MEA  deformation  are  also  affected.  Since  the 
relative  position  between  the  top  and  bottom  end  plate  of  PEM  fuel 
cell  stack  is  fixed  after  assembly,  the  polymer  membrane  is  spa¬ 
tially  confined  under  the  BPP  (with  the  gas  flow  channels)  and  the 
porous  carbon  electrodes,  as  shown  in  Fig.  1.  As  the  RH  increases, 
membrane  absorbs  water,  swells  and  pushes  the  electrodes.  As  a 
consequence,  GDL  will  be  further  compressed  under  the  land  and 
the  protrusion  into  channel  increases  due  to  the  tendency  of  mem¬ 
brane  swelling.  This  membrane  swelling  also  changes  the  local 
contact  forces  due  to  the  redistribution  of  the  stress  field  in  fuel 
cell  stacks. 
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Fig.  1.  Schematic  of  assembly  and  swelling  in  a  PEM  fuel  cell. 


Table  1 

Simulation  cases. 


Table  4 

Young’s  modulus  (MPa)  at  various  temperature  and  humidity  for  Nation®  112  [6]. 


A  sequential  approach  is  implemented  in  this  study.  A  finite 
element  model  (FEM)  is  first  developed  to  model  the  stack  defor¬ 
mation  under  different  levels  of  assembly  pressure,  temperature 
and  RH.  Component  deformation,  the  change  of  material  proper¬ 
ties  and  local  contact  pressure  are  obtained  from  the  FEM  model. 
Then  gas  flow  and  diffusion,  chemical  reactions,  ion  and  electron 
transport  are  modeled  based  the  updated  geometry  and  material 
properties.  Contact  resistance  is  also  analyzed  using  the  model.  The 
impact  of  assembly  pressure  and  operating  conditions  is  evaluated 
by  fuel  cell  performance  and  current  density  distributions. 


Base  case 

Case  1 

Case  2  Case  3 

Pressure  (MPa)  0 

3 

3  3 

T(°C)  25 

25 

85  85 

RH (%)  40 

40 

40  90 

Table  2 

Geometric  and  physical  parameters  for  the  structural  deformation  and  mass  transfer 
analysis. 

Parameter 

Value 

BPP  thickness  (hi ) 

2  mm 

GDL  thickness  (h2) 

100  |jim 

Catalyst  layer  thickness  (h3) 

20  |xm 

Membrane  thickness  (/14) 

50  |jim 

Channel  height  (h5) 

0.5  mm 

Channel  width  (wi ) 

1  mm 

Land  width  (w2) 

1  mm 

Channel  length 

5  mm 

GDL  initial  porosity 

0.6  [14] 

Catalyst  layer  porosity 

0.06  [14] 

GDL  electric  conductivity 

Uncompressed 

In-plane 

3.4  x  104  S  m-1  [12] 

Through-plane 

1.4  x  102  Sm_1  [12] 

GDL  permeability 

1.76  x  10-11  m2 

Catalyst  layer  electronic  conductivity 

lOOSm-1  [15] 

Catalyst  layer  ionic  conductivity 

1.7Sm_1  [16] 

This  paper  investigates  the  influence  of  operating  conditions 
(temperature  and  RH)  on  stack  deformation  and  contact  resistance 
by  improving  the  model  previously  developed  by  the  authors  [2] 
with  respect  to  the  effects  of  assembly  pressure  on  fuel  cell  per¬ 
formance  by  incorporating  temperature  and  RH  effects  with  GDL 
deformation  and  contact  resistance.  Furthermore,  this  paper  also 
investigates  the  current  density  distribution.  During  the  operation 
of  a  PEM  fuel  cell,  significant  variation  of  the  local  current  density 
exists  across  the  plane  of  the  cell.  This  causes  sharp  local  tem¬ 
perature  and  stress  gradient  as  well  as  degrading  the  efficiency  of 
water  management  [8-10].  Current  density  distribution  is  also  an 
important  measure  to  evaluate  fuel  cell  performance  and  durabil¬ 
ity.  Current  density  distributions  under  various  assembly  pressure 
and  operating  conditions  are  investigated. 


2.  Model  description 

A  FEM  based  structural  model  is  developed  to  simulate  stack 
deformation  under  various  assembly  pressures,  temperatures  and 
RHs.  Upon  obtaining  the  deformed  geometry  and  material  proper¬ 
ties  of  GDL  and  membrane,  a  computational  fluid  dynamic  (CFD) 
based  fuel  cell  performance  model  is  developed  to  analyze  the 
multi-component  gas  transport,  chemical  reactions,  charge  transfer 
and  contact  resistance  based  on  the  deformed  GDL  shape  and  mod¬ 
ified  GDL  gas  transport  parameters.  Specifically,  the  local  contact 
force  between  BPP  and  GDL  can  be  obtained.  The  contact  resis¬ 
tance  is  then  simulated  based  on  Zhou  et  al.  [  11  ]  and  included  in  the 
multi-physics  performance  model  to  predict  fuel  cell  performance. 

Four  different  cases  are  modeled  to  analyze  assembly  pressure, 
temperature  and  relative  humidity  impacts  as  shown  in  Table  1 . 

2  A.  Stack  deformation  model  under  elevated  temperature  and 
humidity 

The  model  used  in  the  current  investigation  is  an  extension 
of  the  model  developed  by  Zhou  et  al.  [2].  In  the  current  work, 
temperature-  and  humidity-dependent  properties  of  the  mem¬ 
brane  are  incorporated  in  investigating  stack  deformation  under 
various  assembly  pressures  and  operating  conditions.  The  geomet¬ 
ric  parameters  and  physical  properties  of  the  components  are  listed 
in  Tables  2-4,  where  the  elastic  constants,  coefficients  of  swelling 
and  thermal  expansion  of  the  component  materials,  are  collected 
from  various  references  [6,12,13]. 

A  simple  while  representative  2D,  single-channel,  half-cell  FEM 
model  for  the  cathode  side  is  built  for  stack  deformation  analysis. 
Such  a  simplification  is  pursued  by  taking  advantage  of  the  geomet¬ 
rical  periodicity  of  the  cell,  and  considering  the  fact  that  the  length 
of  gas  channels  is  typically  much  larger  (by  ~100  times)  than  their 
cross-section  dimensions,  which  justifies  the  plane-strain  assump¬ 
tion. 

The  material  properties  for  the  graphite  plates  are  set  to  those 
of  commercial  graphite  material  and  for  the  carbon  paper  from 
TORAY®  TGP-H-030.  It  is  assumed  that  they  do  not  swell  in  response 
to  moisture.  Especially  for  GDL,  it  is  treated  as  nonlinear  elastic 


Table  3 

Component  material  mechanical  properties  [6,12,13]. 


546 


Y.  Zhou  et  al.  /  Journal  of  Power  Sources  192  (2009)  544-551 


Relative  Humidity  (%) 


Fig.  2.  Swelling  expansion  as  a  function  of  humidity  and  temperature  for  Nafion® 
112  [6]. 


[12]  since  its  material  properties  will  settle  down  after  the  first 
few  compression  cycles.  The  compression  curve  for  GDL  is  dif¬ 
ferent  from  decompression  curve.  The  compression  strain-stress 
curve  has  been  applied  to  calculate  GDL  deformation  and  contact 
resistance  since  in  the  model  the  assembly  pressure  on  the  stack  is 
gradually  increased  in  this  analysis. 

Linear  elastic,  perfectly  plastic  constitutive  behavior  with 
temperature-  and  humidity-dependent  material  properties  is 
assumed  for  the  membrane.  Young’s  modulus  and  yield  strength 
of  the  membrane  are  defined  for  four  temperature  and  RH  values 
based  on  experimental  data  for  Nafion®  11 2  from  two  literature. 
Even  though  a  slight  anisotropy  in  the  material  properties  was 
observed  experimentally  [6],  in  this  study  we  assume  for  simplicity 
that  the  material  properties  are  isotropic  (Fig.  2). 

The  model  is  built  using  the  commercial  FEM  software  ABAQUS. 
Four-node  quadrilateral  plane-strain  elements  (CPE4)  are  used  to 
mesh  the  components,  as  shown  in  Fig.  3.  The  bottom  membrane 
surface  is  fixed  vertically  (in  Y-direction),  and  X-symmetry  condi¬ 
tions  are  applied  to  the  side  vertical  boundaries  of  all  components. 
The  component  interfaces  are  bonded,  with  no-slip  allowed.  The 
analysis  consists  of  two  sequential  steps,  as  described  in  the  fol¬ 
lowing. 


2.1.1.  Simulation  of  assembly -induced  stack  deformation 

In  the  stack  assembly,  it  is  desirable  to  have  a  uniform  assembly 
pressure  in  the  BPP,  even  though  the  clamping  load  may  be  local¬ 
ized.  As  such,  a  uniform  assembly  pressure  is  assumed  for  the  BPP. 
Under  ambient  conditions,  i.e.,  25  °C  temperature  and  40%  RH,  a 
series  of  assembly  pressures  of  0.3, 0.5, 0.7, 1.0, 2.0, 3.0,  and  6.0  MPa 
are  applied  to  the  BPP  top  surface.  The  stack  deformation  is  calcu¬ 
lated  under  these  assembly  pressures. 


2.1.2.  Simulation  of  stack  deformation  due  to  swelling  and 
thermal  expansion 

After  the  equilibrium  state  of  the  components  is  achieved  under 
the  applied  pressure,  the  top  surface  of  the  BPP  is  fixed  in  the  ver¬ 
tical  direction  at  the  deformed  position.  Temperature  and  RH  of 
the  whole  stack  are  then  raised  to  the  desired  levels,  i.e.,  85  °C  and 
40%  RH,  and  90%  RH,  under  which  the  components  expand  and 
swell.  Meanwhile,  the  increase  of  temperature  and  RH  also  changes 
the  material  properties,  and  leads  to  the  redistribution  of  deforma¬ 
tions  in  the  whole  stack  under  the  constraint  of  end  plates,  or  the 
constraint  of  BPP  for  the  single-cell. 

The  total  strain  in  the  membrane  due  to  moisture  and  temper¬ 
ature  change  is  calculated  based  on  the  coefficients  of  swelling 
and  thermal  expansion.  The  swelling  behavior  is  temperature- 
dependent,  which  makes  the  total  strain  calculation  difficult  to 
implement  in  FEM  modeling  since  ABAQUS,  which  is  similar  to 
most  of  the  commercial  software  package,  can  only  simulate  the 
expansion  caused  by  the  temperature  field  change.  Thus,  a  new 
parameter,  denoted  as  equivalent  coefficient  of  expansion,  a',  is 
defined  to  overcome  this  challenge  by  combining  the  effects  of 
swelling  and  thermal  expansion.  As  the  temperature  and  humid¬ 
ity  are  increased  to  T  and  RH,  respectively,  a'  can  be  expressed 
as, 


a'(T,  RH)  = 


[1  +  c*(T)AT][l  +  P( RH,  T)ARH(T)]  -  1 
AT 


(1) 


where  a  is  the  coefficient  of  thermal  expansion  (CTE),  ft  is 
the  coefficient  of  swelling  expansion  and  calculated  by  data 
obtained  from  Fig.  2.  ARH  and  AT  are  the  changes  of  RH  and 
temperature,  respectively,  from  one  state  to  another.  When  the 
humidity  effect  is  neglected,  i.e.,  /3(RH,T)  =  0,  a'(T,RH)  =  a(T)  and 
when  the  thermal  expansion  effect  is  neglected,  i.e.,  a(T,RH)  =  0, 
a'(T,RH)  =  /3(RH,T)ARH(T)/AT. 

The  elevated  temperature  and  humidity  cases  simulated  are 
85  °C  and  40%  RH,  and  85  °C,  90%  RH,  which  will  be  compared 
with  the  base  case  at  25  °C  and  40%  RH.  Then  the  change  of  RH 
is  converted  to  the  equivalent  of  temperature  change  using  Eq.  (1 ). 
Through  this  way,  the  swelling  and  thermal  expansion  can  be  both 
considered  in  ABAQUS  simulation. 

Under  the  applied  assembly  pressure  and  the  given  temperature 
and  RH,  the  volumetric  strain  of  every  element  or  the  thickness  of 
the  GDL  can  be  used  to  estimate  its  modified  porosity.  Based  on 
the  deformed  configuration  and  the  modified  GDL  porosity,  mass 
transfer  analysis  is  conducted. 

It  is  assumed  that  the  change  in  thickness  under  compression 
is  due  to  change  in  volume  of  void  space,  not  in  volume  of  solid 
material.  Thus,  the  porosity  of  the  GDL  can  be  calculated 


0/  „  v'  ^solid 

£  =  £  o - 

^0  —  ^solid 


(2) 


where  £0  is  the  initial  porosity  of  GDL,  V0  is  the  uncompressed  vol¬ 
ume  and  V  is  the  volume  after  compression.  From  the  structure 
model  above,  the  individual  element  volume  can  be  exported  so 
that  local  porosity  is  calculated  accordingly. 
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Assembly  Pressure  (MPa) 

Fig.  4.  Contact  resistance  vs.  assembly  pressure. 

2.2.  Contact  resistance 

Contact  resistance  constitutes  a  significant  part  of  ohmic  resis¬ 
tance  in  the  fuel  cell,  and  needs  to  be  considered  when  evaluating 
fuel  cell  performance.  A  high  contact  resistance  exists  between  the 
GDL  and  BPP  due  to  the  fact  that  the  GDL  is  porous  (60%  to  80% 
porosity)  and  the  surface  of  BPP  engaging  GDL  is  rough.  Hence,  cur¬ 
rent  flow  between  BPP  and  GDL  occurs  only  at  sites  where  the  BPP 
contacts  GDL. 

Assembly  pressure  affects  contact  resistance  in  PEM  fuel  cells. 
Contact  resistance  can  be  estimated  based  on  surface  roughness 
parameters  of  BPP  and  features  of  GDL  structure.  A  detailed  descrip¬ 
tion  of  the  model  is  presented  in  Zhou  et  al.  [11  ]. 

BPP  surface  topology  is  simulated  as  randomly  distributed 
asperities,  and  is  based  on  measured  surface  roughness.  The  GDL  is 
modeled  as  randomly  distributed  cylindrical  fibers.  Upon  obtaining 
these  two  simulated  surfaces,  each  contact  spot  is  located  according 
to  relative  positions.  The  total  resistance  and  pressure  are  obtained 
by  considering  all  contact  spots  as  resistances  in  parallel  and  sum¬ 
ming  the  results  together. 

In  this  study,  the  contact  resistance  between  graphite  BPP  and 
carbon  fiber  paper  GDL,  the  most  common  materials  used  in  PEM 
fuel  cells,  is  applied  in  the  model.  The  surface  roughness  parameters 
obtained  from  the  average  values  of  several  scans  are:  peak  den¬ 
sity  Dpeak=150mm_1,  mean  asperity  summit  radius  Ri=3.26fxm, 
and  variance  of  the  summit  height  distribution  ors=0.728  p,m.  The 
model  predicted  results  of  contact  resistance  change  with  assembly 
pressure  is  shown  in  Fig.  4. 

Since  the  operating  conditions  will  influence  local  contact  force 
and  GDL  porosity,  two  crucial  factors  for  determining  contact  resis¬ 
tance,  the  changes  of  contact  resistance  should  be  considered.  The 
contact  force  and  GDL  porosity  will  be  obtained  from  the  structural 
model  in  Section  2.1. 

2.3.  PEM  fuel  cell  performance  analysis 

A  FEM  computational  fluid  dynamics  package,  COMSOL  Multi¬ 
physics®,  is  applied  to  solve  the  coupled  nonlinear  equations 
representing  the  physical  phenomena  of  gas  and  charge  transfer. 
Fig.  5  shows  the  model  domain  for  the  base  case  with  no  assem¬ 
bly  pressure,  25  °C  temperature  and  40%  RH.  The  3D  model  domain 
consists  of  a  conventional  gas  channel,  cathode  GDL,  catalyst  layer 
(CL),  membrane,  and  anode  electrode  interface.  Because  of  the 
symmetrical  characteristic  of  the  structure,  only  half  the  width 


Fig.  5.  Computational  domain  for  the  base  case. 


of  the  shoulder  and  channel  is  analyzed  to  reduce  computational 
time. 

There  exist  coupled  gas  mass  transfer  and  charge  transfer  in 
PEM  fuel  cell  stacks.  In  the  BPP  channel,  multi-component  gases 
consist  of  N2,  02  and  H20.  Gas  flow  and  components  in  BPP  chan¬ 
nel  is  dependant  on  diffusion  in  GDL  and  oxygen  consumption 
rate  at  the  catalyst  layer.  The  gas  flow  is  modeled  in  combination 
with  diffusion  and  convective  transport.  In  the  GDL  and  CL,  N2,  02 
and  H20  transport  is  driven  by  diffusion  mechanism.  Three  sets 
of  Navier-Stroke  equations  are  used  to  model  the  incompressible 
gas  flow  of  N2,  02  and  H20  in  channels  and  another  three  sets  of 
Maxwell-Stefan  equations  are  used  to  simulated  multi-component 
diffusion  and  convection  in  channels  and  gas  distribution  layers. 
The  velocity  of  gases  through  GDL  is  modeled  as  Darcy’s  flow  since 
there  exists  pressure  difference  across  the  GDL.  The  rate  of  02  con¬ 
sumption  and  H20  production  are  determined  by  coupled  effects 
of  gas  transfer  rate  and  charge  transfer  processes.  Those  coupled 
equations  will  be  solved  simultaneously  to  obtain  the  steady  state 
of  gas  flow,  diffusion  and  chemical  reactions. 

Charge  balance  in  the  solid  phase  (GDL,  membrane  and  CL)  is 
taken  into  account  in  this  model,  including  the  ion  transport  in  CL 
and  membrane  and  electron  transport  in  GDL  and  CL.  In  GDL,  due  to 
no  chemical  reactions,  the  electronic  charge  is  conserved.  Within 
the  CL,  electrons  are  consumed  by  the  02  reduction  reaction  and 
the  02  consumption  rate  varies  across  the  CL.  The  gradient  of  the 
electronic  current  is  proportional  to  the  oxygen  consumption  rate. 
The  contact  resistance  between  GDL  and  the  land  of  BPP  is  taken 
into  account  by  modeling  the  boundary  as  a  thin  layer.  The  electric 
conductivity  and  thickness  of  the  boundary  layer  are  calculated  so 
that  its  total  resistance  corresponded  to  the  contact  resistance  val¬ 
ues  obtained  from  the  mechanical  deformation  model  and  contact 
resistance  simulation  in  Section  2.2.  Ion  transfer  only  takes  place 
in  CL  and  membrane.  Ions  are  consumed  at  the  CL  due  to  chemical 
reactions.  The  ion  transfer  is  conserved  in  the  membrane  domain. 

For  boundary  conditions,  gas  mixture  is  assumed  to  enter  the 
channel  normal  to  the  inlet  cross-section.  All  walls  in  the  channel 
have  no-slip  boundary  conditions.  The  mass  and  momentum  trans¬ 
port  boundary  conditions  between  the  BPP  shoulders  and  the  GDL 
are  all  insulated.  At  the  interface  between  GDL  and  CL,  it  is  assumed 
that  no  contact  resistance  exists.  The  electronic  current  and  the 
fluxes  of  N2, 02  and  H20  in  the  Y-direction  are  continuous.  The  flux 
of  protons  is  set  to  zero  because  there  is  no  ionic  phase  in  the  GDL. 
At  the  boundary  between  the  CL  and  the  membrane,  the  fluxes  of 
N2,  02  and  H20,  and  electronic  current  in  the  Y-direction  are  set  to 
zero  while  the  ionic  current  are  continuous.  The  anode  is  assumed 
to  be  an  interface,  and  the  ionic  potential  is  approximately  zero  at 
the  interface  between  the  anode  CL  and  the  membrane  based  on 
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Fig.  6.  Model  predicted  contour  of  displacement  at  Y-direction  under  3  MPa  assembly  pressure. 


the  assumption  that  the  hydrogen  oxidation  reaction  rate  is  so  fast 
that  the  anodic  overpotential  is  negligible  [14].  All  other  boundaries 
are  assumed  to  be  insulated  due  to  symmetry. 

Consistent  with  traditional  fuel  cell  models  [17,18],  the  model 
is  assumed  to  be  steady  state,  isothermal  and  isobaric.  Another 
assumption  is  that  single-phase  water  transport.  Water  exists  only 
in  the  gas  phase.  The  existence  of  liquid  water  is  taken  into  account 
by  using  effective  porosity  as  the  initial  porosity  for  GDL,  which  may 
be  significantly  smaller  than  the  raw  material  because  of  possible 
flooding  by  liquid  water. 

The  geometric  and  physical  properties  of  the  components  are 
listed  in  Table  2.  The  electrical  conductivity  in  the  CL  is  assumed  to 
be  lOOSm-1  for  all  the  cases.  For  the  GDL,  the  electrical  conduc¬ 
tivity  is  anisotropic.  The  relationship  between  in-plane  oxz  vs.  the 
through-plane  conductivity  oy  follows  the  relationship  related  to 
channel  geometry  and  thickness  of  GDL  [12]: 

Oxz  =  (Wi  +W2)Wi 

Oy  8  b\ 

A  wide  range  of  GDL  conductivity  values  have  been  employed  in 
various  PEM  fuel  cell  models.  In  this  study,  based  on  the  experi¬ 
mental  data  from  Mathias  et  al.  [12],  through-plane  conductivity 
for  Toray  carbon  fiber  paper  is  1369  S  nrr 1  in  the  base  case  (no  com¬ 
pression).  The  in-plane  conductivity  of  carbon  fiber  paper  is  over 
an  order  of  magnitude  larger  than  the  through-plane  value.  And 
the  channel/land  geometry  imposes  stricter  requirements  on  the 
through-plane  conductivity. 

Both  the  in-plane  and  through-plane  conductivities  increase 
monotonically  with  increase  of  assembly  pressure,  with  the  change 
in  through-plane  conductivity  especially  significant.  The  change  of 
GDL  conductivities  follows  the  relationship  obtained  experimen¬ 
tally  by  [4],  as  shown  in  Table  5. 

Some  parameters  change  with  temperature  and  humidity.  The 
open  circuit  potential  E°  for  the  overall  reaction  is  calculated  as 
[19]: 

E°  =  0.2329  +  0.00251  (4) 

The  binary  diffusivities  Dijt  obtained  experimentally  at  atmo¬ 
spheric  pressure,  are  scaled  with  the  temperature  and  pressure 
according  to  Cussler  [20]. 

D».Dli(7i.p„)iQ'-5  (5) 

The  membrane  conductivity  depends  on  RH.  The  following  rela¬ 
tionship  between  ionic  conductivity  a  (in  the  unit  of  Scnrr1)  and 
RH  obtained  from  curve  fitting  the  experimental  data  is  employed 
[21]. 

a  =  1.3  x  10“7exp(14RHa2)  (6) 


Table  5 

GDL  electrical  conductivities. 


GDL  electrical  conductivity 

Base  case  (uncompressed) 

3  MPa 

In-plane 

3.4  x  104  S  m~ 

1  [121 

3.43  x  104  S  m_1 

Through-plane 

1.4  x  102  S  m-1 

[12] 

l.eSxKPSrrr1 

Using  deformed  geometry  and  parameters  for  each  case,  the  rela¬ 
tionship  between  cell  voltage  and  current  density  can  be  obtained 
by  setting  the  potential  of  cathode  current  collector  to  cell  volt¬ 
age.  COMSOL  Multi-physics®  is  used  to  solve  the  coupled  nonlinear 
equations  for  gas  and  charge  transfer.  Details  of  the  model,  geom¬ 
etry  meshing  and  solver  settings  were  consistent  with  previous 
model  [1]. 

3.  Results  and  discussion 

3.1.  GDL  and  membrane  deformation 

The  membrane  and  GDL  are  pre-stressed  due  to  the  applied 
clamping  as  explained  previously  The  displacement  under  3  MPa 
of  assembly  pressure  (25  °C  temperature  and  40%  RH)  is  calculated 
and  shown  in  Fig.  6.  Compared  to  GDL  and  membrane,  BPP  is  a  very 
stiff  material  and  has  very  small  deformation.  Fig.  6  only  shows  the 
bottom  parts  of  Fig.  3,  including  GDL  and  membrane.  The  GDL  is 
deformed  severely  under  the  land,  but  it  is  almost  unchanged  under 
the  channel.  The  membrane  has  very  little  deformation  because  it 
has  larger  Young’s  modulus  compared  to  that  of  GDL. 

Fig.  7  shows  the  deformation  of  GDL  and  membrane  when  tem¬ 
perature  and  RH  change  under  the  same  assembly  pressure  (3  MPa). 
When  temperature  and  RH  increase,  swelling  of  membrane  pushes 
GDL  more  into  the  channel  cavity  and  the  contact  surface  between 
GDL  and  membrane  becomes  curved.  But  the  thickness  of  GDL 
under  the  channel  almost  remains  the  same  and  that  the  total 
change  from  the  original  uncompressed  volume  remains  small. 

BPP,  GDL  and  membrane  all  have  thermal  expansion  to  some 
extent,  although  the  thermal  expansion  coefficients  for  GDL  and 
BPP  are  very  small.  For  membrane,  the  special  characteristic  is  that 
it  also  expands  with  respect  to  humidity.  Therefore,  the  total  strain 
in  the  membrane  due  to  change  in  moisture  and  temperature  is 
calculated  from  the  coefficients  of  swelling  and  thermal  expansion. 
When  temperature  increases  to  85  °C,  there  is  a  slight  variation  of 


X  (mm) 

Fig.  7.  Deformed  shape  of  GDL  and  membrane  at  various  conditions. 
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3MPa/25°c  /40%  3MPa/85°C  /40%  3MPa/85°C  /90% 


Fig.  8.  Contact  pressure  and  contact  resistance  at  various  operation  conditions. 

the  GDL  and  membrane  deformation.  It  is  mainly  due  to  the  fact 
that  the  temperature  induced  strains  for  both  GDL  and  membrane 
are  very  small  even  though  the  Young’s  modulus  of  membrane 
drops  by  about  40%.  Swelling  expansion  coefficient  due  to  mois¬ 
ture  absorption  is  defined  as  the  relative  change  in  length  per  1% 
RH  change.  When  RH  increases  to  90%,  membrane  becomes  softer 
since  its  Young’s  modulus  drops  from  197  MPa  to  46  MPa.  Whereas, 
swelling  induced  strain  for  membrane  is  significant. 

Fig.  8  shows  the  change  of  contact  pressure  and  contact  resis¬ 
tance.  The  contact  pressure  at  3  MPa/25  °C/40%RH  condition  is 
5.78  MPa  and  contact  resistance  is  0.22  mC2  cm2.  When  the  temper¬ 
ature  increases  to  85  °C,  the  contact  pressure  goes  up  to  7.2  MPa  and 
contact  resistance  decreases  to  0.198  mC2  cm2.  When  RH  further 
increases  from  40%  to  90%  at  85  °C,  the  contact  pressure  increases 
to  9.3  MPa  and  contact  resistance  is  0.12  mC2  cm2.  That  is  due  to  the 
re  distribution  of  the  stress  field  in  the  stack.  Since  the  relative  loca¬ 
tion  between  upper  surface  of  BPP  and  lower  surface  of  membrane 
is  fixed  after  assembly,  the  resulting  displacement  at  the  upper 
boundary  from  the  previous  step  is  fixed  throughout  the  analysis. 
The  internal  forces  have  to  be  redistributed  to  accommodate  addi¬ 
tional  loading.  Hence  the  contact  force  is  changed  when  fuel  cell 
operating  conditions  change.  When  membrane  swells,  the  contact 
force  becomes  larger.  The  contact  force  between  BPP  shoulder  and 
GDL  determines  the  contact  resistance  in  between.  The  change  of 
contact  resistance  should  be  considered  when  evaluating  fuel  cell 
overall  performance. 

When  the  stack  is  not  well  assembled,  the  contact  resistance 
could  be  very  high,  ranging  up  to  a  few  hundreds  mC2  cm2  in  a  fuel 
cell  stack  [22].  A  reasonable  value  of  20  cm2  is  assumed  in  the 
base  case  model  [23].  The  contact  resistance  under  3  MPa  assembly 
pressure  and  various  operating  conditions  are  significantly  smaller 
due  to  higher  electrical  conductivity  and  smoother  contacting  sur¬ 
faces  used  in  the  simulation. 

3.2.  Performance  model 
3.2.1.  Polarization  curves 

The  polarization  curves  are  obtained  by  solving  the  average 
current  density  for  different  values  of  cell  voltage.  Fig.  9  shows 
the  polarization  curves  for  the  four  cases.  In  brief,  assembly  pres¬ 
sure  causes  performance  decrease  but  higher  temperature  and  RH 
increase  the  polarization  curves  over  the  entire  operating  range  of 
the  cell. 


It  is  observed  that  under  the  assembly  pressure  of  3  MPa,  the 
performance  decreases  which  is  consistent  with  previous  findings 
[1  ].  Such  a  performance  deterioration  is  mainly  due  to  the  fact  that 
compression  causes  a  decrease  of  GDL  porosity  which  then  imposes 
more  resistance  to  gas  flow,  even  though  compression  reduces  con¬ 
tact  resistance.  There  are  no  significant  differences  between  the 
compression  case  (3  MPa/25  °C/40%RH)  and  the  base  case  at  the 
cell  voltage  level,  indicating  that  the  overall  cell  performance  is  not 
significantly  affected  by  the  assembly  pressure  at  3  MPa  for  this  set 
of  component  geometry  and  material  properties. 

When  temperature  increases  from  25  °C  to  85  °C,  the  perfor¬ 
mance  curve  has  a  moderate  increase.  This  is  due  to  the  total  effects 
of  contact  resistance  reduction,  GDL  protrusion  and  the  change  of 
some  kinetic  parameters.  Contact  resistance  reduction  leads  to  the 
decrease  of  ohmic  voltage  loss.  GDL  deformation  reduces  gas  flow 
area  and  through-plane  thickness,  which  can  facilitate  the  gas  flow 
therefore  reduces  flow  resistance.  Also,  the  diffusion  characteristics 
of  gas  mixture  and  kinetic  parameters  increase  with  temperature, 
which  benefit  for  gas  diffusion,  and  result  in  the  decrease  of  the 
cathode  activation  overpotential  with  temperature. 

When  RH  goes  to  90%,  the  performance  curve  has  a  more  sig¬ 
nificant  increase.  When  RH  increases,  the  ionic  conductivity  the 
membrane  increases,  so  that  the  overpotential  is  decreased.  In 
terms  of  cell  performance,  this  is  desirable  because  lower  over¬ 
potentials  result  in  higher  cell  voltage.  Also  the  change  of  RH  has 
more  impacts  on  contact  pressure  and  contact  resistance  between 
BPP  and  GDL,  as  shown  in  Fig.  8.  When  RH  increases  from  40%  to 
90%,  the  contact  resistance  drops  more.  This  is  one  of  the  contribut¬ 
ing  factors  for  the  performance  increase.  In  addition,  GDL  is  more 
protruded  into  BPP  channel  due  to  swelling  and  the  change  of  mate¬ 
rial  properties,  and  the  change  of  porosity  and  flow  channel  area 
reduction  are  beneficial  to  gas  flow  and  diffusion,  which  will  in  turn 
improve  performance.  All  those  effects  integrate  together  to  yield 
the  performance  better  with  higher  temperature  and  humidity. 

3.2.2.  Current  production  rate  in  CL  along  the  channel  length 
direction 

Fig.  10(a)  and  (b)  presents  the  current  production  rate  and  vari¬ 
ation  at  the  interface  between  GDL  and  catalyst  layer  along  the 
channel  direction  for  three  cases  when  the  cell  voltage  is  0.5  V. 
When  the  assembly  pressure,  elevated  temperature  and  RH  are 
all  taken  into  account,  there  is  a  significant  change  of  the  current 
production  rate  at  this  interface. 
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Fig.  10.  Current  production  rate  profiles  at  the  interface  between  CL  and  GDL  for  different  cases  at  0.5  V. 
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When  an  assembly  pressure  of  3  MPa  is  applied,  the  average 
current  density  of  the  case  at  3  MPa/25  °C/40%  is  lower  than  the 
base  case  (0  MPa/25  °C/40%)  since  assembly  pressure  causes  perfor¬ 
mance  decrease.  The  current  production  rate  starts  to  have  some 
variations  along  the  channel  length  direction  at  this  interface.  At 
the  base  case,  the  current  production  rate  only  has  1.7%  varia¬ 


tion  compared  with  the  average  value  over  the  whole  surface. 
At  the  compression  case  (3  MPa/25  °C/40%),  the  current  produc¬ 
tion  rate  variation  is  increased  to  10%.  This  is  due  to  the  fact 
that  assembly  pressure  induces  large  unevenness  of  GDL  poros¬ 
ity.  GDL  porosity  is  approximately  0.3  under  the  BPP  land  while 
0.6  under  the  BPP  channel  for  the  compression  case.  This  vari- 
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Fig.  11.  Current  production  rate  profiles  in  the  CL  for  different  cases  at  0.5  V. 
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ation  in  GDL  porosity  causes  the  variation  in  current  production 
rate. 

When  temperature  rises  to  85  °C,  the  average  current  production 
rate  is  higher  since  the  voltage  output  increases  with  temperature, 
while  the  variation  of  current  production  rate  increases  to  15%.  Tem¬ 
perature  increase  causes  more  inhomogeneous  compression,  which 
in  turn  causes  more  current  generation  under  the  channel  and  less 
current  generation  under  the  land. 

When  RH  increases  to  90%,  the  average  current  production  rate  is 
substantially  increased  since  the  voltage  output  increases  greatly  at 
higher  RH  because  of  higher  membrane  conductivity  and  smaller 
contact  resistance.  Meanwhile,  the  current  production  rate  vari¬ 
ation  is  also  increased  to  68%.  A  notable  portion  of  current  is 
generated  mostly  under  the  channel  and  the  current  production 
rate  is  reduced  under  the  BPP  land.  That  is  attributed  by  the  fact 
that  the  inhomogeneous  compression  is  much  more  severe  when 
RH  increases  as  shown  in  Fig.  7.  Thus  GDL  porosity  unevenness 
becomes  more  significant.  It  is  more  difficult  for  the  gases  to  trans¬ 
port  into  the  GDL  under  the  land;  hence  the  GDL  under  the  land  is 
in  the  starvation  of  reactant  gases  and  generates  less  current. 

3.2.3.  Current  production  rate  in  CL  along  the  thickness  direction 

Fig.  11  shows  the  distributions  of  current  production  rate  in  the 
thickness  direction  of  the  CL  at  the  cell  voltage  of  0.5  V  and  Z  =  2  mm. 
The  same  trends  as  in  Fig.  10  can  be  observed.  Both  compression  and 
elevated  temperature  and  RH  cause  significant  variation  along  the 
CL  thickness  and  width  direction.  Less  current  is  produced  under 
the  land  area,  because  of  the  reduced  oxygen  diffusion  in  com¬ 
pressed  region.  More  current  is  generated  under  the  channel  area 
which  is  caused  by  the  combined  effects  including  the  dependence 
of  diffusion  parameters  and  kinetic  parameters  on  temperature  and 
RH,  GDL  protrusion  and  porosity  variation.  As  shown  in  Fig.  7,  GDL  is 
expanded  under  the  channel  area  so  that  the  porosity  is  increased. 
Thus  the  gas  transfer  is  facilitated  and  more  current  is  generated. 

4.  Conclusions 

In  this  paper  we  presented  a  comprehensive  analysis  of  the  influ¬ 
ence  of  elevated  temperature  and  relative  humidity  on  fuel  cell 
stack  deformation,  contact  resistance  and  performance.  Elevated 


temperature  and  humidity,  especially  humidity,  caused  severe  GDL 
protrusion  and  stress-redistribution,  which  lead  to  the  reduction 
of  contact  resistance.  Even  though  elevated  temperature  and  RH 
improved  fuel  cell  performance,  significant  variation  of  current  pro¬ 
duction  rate  was  observed  along  all  directions  in  the  cell  when 
assembly  pressure  was  applied  and  temperature  and  RH  were 
increased  close  to  the  operation  conditions.  The  stack  would  be 
more  prone  to  degradation  with  such  a  significant  variation.  One  of 
the  future  research  topics  is  the  assembly  optimization  at  various 
operating  conditions  to  improve  performance  as  well  as  mitigate 
the  current  variation. 
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